Simulations of Binary Black Hole Mergers Using Spectral Methods 
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Several improvements in numerical methods and gauge choice are presented that make it possible 
now to perform simulations of the merger and ringdown phases of "generic" binary black-hole 
evolutions using the pseudo-spectral evolution code SpEC. These improvements include the use of 
a new damped-wave gauge condition, a new grid structure with appropriate filtering that improves 
stability, and better adaptivity in conforming the grid structures to the shapes and sizes of the 
black holes. Simulations illustrating the success of these new methods are presented for a variety of 
binary black-hole systems. These include fairly "generic" systems with unequal masses (up to 2:1 
mass ratios), and spins (with magnitudes up to QAM 2 ) pointing in various directions. 



I. INTRODUCTION 



Black-hole science took a great stride forward in 2005 
when Pretorius [l[ performed the first successful full non- 
linear dynamical numerical simulation of the inspiral, 
merger, and ringdown of an orbiting black-hole binary 
system; this initial success then stimulated other groups 
to match this achievement within months 0, H[ . These 
developments lead quickly to advances in our understand- 
ing of black-hole physics: investigations of the orbital 
mechanics of spinning binaries [1, HI, H, 0, @, @ , studies of 
the recoil from the merger of unequal mass binary sys- 
tems E2, El, EH , the remarkable discovery of unex- 
pectedly large recoil velocities from the merg er of certain 
^inningbinarysystems @, EI EE El, E3, El, El, HE 
EL HI, IH, [H, HE , investigations into the mapping 
between the binary black-hole initial conditions (individ- 
ual masses and spins) and the final state of the merged 
black hole [H HI, Hi, HE HH, HI, HE HI, and improve- 
ments in our understanding of the validity of approxi- 
mate binary black-hole orbital calculations u sing post- 
Newtonian methods [H, EE EE EE El, HE HO, Ha ■ 

These first results on binary black-hole systems were 
obtained by several different groups using different codes 
based on two different formulations of the Einstein equa- 
tions (BSSN and generalized harmonic), using two differ- 
ent methods for treating the black-hole interiors (mov- 
ing puncture and excision), and using rather different 
gauge conditions to fix the spacetimc coordinates. All 
of these early results, however, were obtained with codes 
based on finite difference numerical methods and adap- 
tive mesh refinement for computational efficiency. The 
Caltech/Cornell collaboration decided a number of years 
ago to follow a different path by developing an Einstein 
evolution code, called SpEC, based on spectral methods. 
The advantages of spectral methods are their superior ac- 
curacy and computational efficiency, and their extremely 
low numerical dissipation that is ideal for solving wave 
propagation problems with high precision. The disadvan- 
tages are the relative complexity of spectral codes, the 
lack of any appropriate pre-existing spectral code infras- 
tructure (analagous to CACTUS [43[ for example), and 
the extreme sensitivity of spectral algorithms to develop- 



ing instabilities when any aspect of the solution method is 
mathematicaly ill-posed. Consequently it has taken our 
group somewhat longer to bring SpEC up to the level of 
the state-of-the-art codes in the field. 

We, along with our Caltech/Cornell collaborators, 
have developed a large number of numerical and analyti- 
cal tools over the years that make it possible for SpEC to 
evolve binary black-hole systems with great er p recision 
than any other code (at the present time) j44|- These 
technical developments include the derivation and imple- 
mentation of constraint preserving and physical (no in- 
coming gravitational wave) boundary conditions [IE E|| , 
dual frame evolution methods and feedback and control 
systems to lock the computational grids onto the location 
of the black holes [l?) , and special angular filtering meth- 
ods needed to cure an instability that occurs when tensor 
fields are evolved [48|. Using these methods we and our 
collaborators have performed a number of high precision 
evolutions of the inspiral portions of binary black-hole 
systems, and have used the gravitational waveforms from 
these evolutions to calibrate the accuracy of the approxi- 
mate post-Newtonian waveforms that are widely used in 
gravitational wave data analysis [IE EE EH, E3 . 

During the past year our group has begun to have some 
success in performing the more dynamical and difficult 
(for spectral codes) merger and ringdown portions of bi- 
nary black-hole evolutions [H EE M, EI EE HI ■ The 
techniques we used to achieve these first merger calcula- 
tions turned out to be rather non-robust however, requir- 
ing a great deal of hands-on adjustment and fine tuning 
for each new case we attempted. In this paper we report 
on several new numerical and analytical developments 
that allow us now to perform stable and accurate binary 
merger and ringdown simulations of a fairly wide range 
of binary black-hole systems. These new methods ap- 
pear to be quite robust, and no fine tuning is required: 
the same basic method works on each of the cases we 
have attempted. 

The three new technical breakthroughs that allow us 
now to perform successful binary black-hole merger and 
ringdown simulations are described in Sec. [II] These ma- 
jor advances include the development of a new gauge 
condition that works extremely well with the generalized 
harmonic formulation of the Einstein system used in our 
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code. This new gauge chooses spatial coordinates that 
are solutions of a damped wave equation, and chooses 
the time coordinate in a way that limits the growth of 
y^/g/N, here g is the determinant of the spatial metric 
and N is the lapse function. This new gauge condition is 
described in some detail in Sec. Ill Al Another important 
development, reported in Sec. Ill Bl is the construction of 
a new grid structure on which our binary black-hole evo- 
lutions are performed. This new non-overlapping grid 
structure, and a certain type of spectral filtering used 
with the new grid structure, removes a class of numerical 
instabilities that were a limiting factor in our ability to 
perform the highly dynamical portions of merger calcula- 
tions, and allows us to increase resolution near the black 
holes in a more targeted and computationally efficient 
way. The third technical development is a new more effi- 
cient and robust method to conform the structure of the 
grid to the shape and size of the black holes in an auto- 
matic dynamical way. These developments are described 
in Sec. Ill CI and the Appendix,, 

Using these new technical tools we have now performed 
successful merger and the succeeding ringdown simula- 
tions of a fairly wide range of binary black-hole systems. 
We describe in Sec. Mil merger simulations for six reason- 
ably diverse cases. These include two systems in which 
the black holes are non-spinning: one has equal mass 
black holes, the other has holes with a 2:1 mass ratio. 
We also describe successful merger and ringdown calcu- 
lations with two different equal-mass binary systems in 
which the holes have identical intrinsic spins of magni- 
tude 0AM 2 : aligned to the orbital angular momentum in 
one system, and anti-aligned in the other. We have also 
performed successful merger and ringdown simulations 
on two fairly generic binary systems. These two systems 
have black holes with the same mass ratio M^/Mi = 2, 
and the same "randomly oriented" initial spins of magni- 
tudes 0.2M 2 and 0.4A/| respectively. But the initial sep- 
arations of the holes are different in the two cases: one 
starts about 8.5 orbits before merger, and another (used 
to perform careful convergence studies) starts about 1.5 
orbits before merger. The same methods were used to 
perform all of these merger and ringdown simulations, 
and no particular fine tuning was required. 



II. TECHNICAL DEVELOPMENTS 

This section describes the three new technical break- 
throughs that allow us now to perform successful binary 
black-hole merger and ringdown simulations. These ma- 
jor advances include the development of a new gauge 
condition, described in Sec. Ill Ai in which the spatial co- 
ordinates satisfy a damped wave equation and the time 
coordinate is chosen in a way that controls the growth 
of the spatial volume element. Another important devel- 
opment, described in Sec. Ill Bl is a new non-overlapping 
grid structure for our binary black-hole simulations, and 
a type of spectral filtering for these new grid structures 



that improves their accuracy and stability. The third new 
technical development, described in Sec. Ill CI is a more 
efficient and robust method of conforming the structure 
of the grid to the shape and size of the black holes in a 
dynamical and automatic way. 

A. Damped- Wave Gauge 

Harmonic gauge is defined by the condition that each 
coordinate x a satisfies the co-variant scalar wave equa- 
tion: 

V c V c a; a = H a = 0. (1) 

Harmonic coordinates have proven to be extremely use- 
ful for analytical studies of the Einstein equations, but 
have found only limited success in numerical problems 
like simulations of complicated highly dynamical black- 
hole mergers. One reason for some of these difficulties 
is the wealth of "interesting" dynamical solutions to the 
harmonic gauge condition itself, Eq. (TTJ). Since all "physi- 
cal" dynamical fields are expressed in terms of the coordi- 
nates, an ideal gauge condition would limit coordinates to 
those that are simple, straightforward, dependable, and 
non-singular; having "interesting" dynamics of their own 
is not a desirable feature for coordinates. The dynamical 
range available to harmonic coordinates can be reduced 
by adding a damping term to the equation: (57l . l58j 

V c V c x a = ^st c d c x a = ^ s t a , (2) 

where t a is the future directed unit normal to the 
constant-t hypersurfaces. Adding such a damping term 
to the equations for the spatial coordinates x % tends to 
remove extraneous gauge dynamics and drives the coor- 
dinates toward solutions of the co- variant spatial Laplace 
equation on the timescale l//i. Choosing l//i to be com- 
parable to (or smaller than) the characteristic timescale 
of a particular problem should remove any extraneous co- 
ordinate dynamics on timescales shorter than the phys- 
ical timescale. The addition of such a damping term in 
the time-coordinate equation is not appropriate however. 
Such a damped-wave time coordinate is driven toward a 
constant value, and therefore toward a state in which it 
fails to be a useful time coordinate at all. It makes sense 
then to use the damped-wave gauge condition only for 
the spatial coordinates: 

V c V c sc' = H l = ust* = -nsN'/N, (3) 

where N l is the shift, and N is the lapse. The appro- 
priate contra- variant version of this damped-wave gauge 
condition is therefore 

H a = -n s g al N l /N, (4) 

where g a \, is the spatial metric of the constant-t hypersur- 
faces. (Note that we use Latin letters from the beginning 
of the alphabet, a, b, c, . . ., to denote four-dimensional 
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spacctime indices, and letters from the middle of the al- 
phabet, i,j, k, . . ., for three-dimensional spatial indices.) 

While the damped-wave gauge is a poor choice for the 
time coordinate, the idea of imposing a gauge that adds 
dissipation to the gauge dynamics of the time coordinate 
is attractive. To find the appropriate expression for t a H a , 
the component of H a not fixed by Eq. 0$, we note that 
the gauge constraint H a + T a = implies that t a H a is 
given by 



(5) 



where g = det gij is the spatial volume element. In our 
experience, a frequent symptom of the failure of simpler 
gauge conditions in binary black-hole simulations is an 
explosive growth of g in the spacetime region near the 
black-hole horizons. This suggests that a good use of the 
remaining gauge freedom would be to attempt to control 
the growth of the spatial volume element, g. Choosing 
the gauge condition, 
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together with Eq. implies the following evolution 
equation for ^fg/N: 



t a d a log 



N 



log 



N 



N^dkN*, (7) 



whose solutions tend to suppress growth in ^/~g/N . The 
discussion of this gauge condition in Ref. [57| shows that 
it also implies that the lapse N satisfys a damped wave 
equation, with the damping factor fi^. So in this sense, 
the gauge condition on the time coordinate, Eq. (|6|), is 
the natural extension of the spatial-coordinate damped- 
wave gauge condition, Eq. ([4]). 

Combining this new lapse condition, Eq. ([5]), with the 
damped-wave spatial coordinate condition, Eq. gives 
the gauge-source function for our full damped-wave gauge 
condition: 



t a 



VsN^ga 



(8) 



The damping factors > and us > can be chosen 
quite arbitrarily as functions of spacetime coordinates x a , 
or even as functions of the spacetime metric V'ob- The 
gauge source function H a depends only on coordinates 
and the spacetime metric ^ a b in this case, so these gauge 
conditions can be implemented directly in the generalized 
harmonic Einstein system without the need for a gauge 
driver. Previous studies of this condition (developed ini- 
tially as a test of the first-order gauge-driver system (5?j|) 
showed it to be quite useful for evolving single black-hole 
spacetimes. In those tests, which included several differ- 
ent evolutions of maximal-slice Schwarzschild initial data 
with large non-spherical gauge perturbations, the black 



hole always evolved quickly toward a non-singular time- 
independent equilibrium state. 

The solutions to the lapse gauge condition, Eq. ([7]), 
can be thought of as equating \og(^/g/N) to a certain 
weighted time average of dkN k /N . The timescale asso- 
ciated with this time averaging is set by ^l, which de- 
termines for example the rate at which y^/iV is driven 
toward an asymptotic equilibrium state. When /i£ is con- 
stant, it is easy to show that Vd/^ 1S driven exponen- 
tially toward this asymptotic state. In highly dynamical 
spacetimes, however, we find that y^/iV must be driven 
even faster than exponential in order to prevent the for- 
mation of singularities in g. This can be accomplished by 
making /i^ larger whenever g becomes large. In practice 
we find that the choice fiL = A*o[log(Y / g/A r )] 2 , (where /zo 
is a constant or perhaps a function of time) is very effec- 
tive at suppressing the growth of these singularities. We 
find that choosing the same damping factor /xg, 
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(9) 



for the spatial part of the gauge condition is also quite 
effective. The binary black-hole merger and ringdown 
simulations described in Sec. IIIII use these and /zg, 
with fiQ taken to be an order-unity function of time (to 
accomodate starting up evolutions from initial data sat- 
isfmg a different gauge condition). 



B. Grid Structure and Filtering 

The pseudo-spectral numerical methods used by our 
code, the Spectral Einstein Code (SpEC) [H, H3, HH , rep- 
resent dynamical fields on a spatial grid structure that is 
specially constructed for each problem. For binary black- 
hole simulations, our group has been using grid structures 
constructed from layers of spherical-shell subdomains 
centered on each black hole, surrounded and connected 
by cylindrical, cylindrical-shell, and/or rectangular-block 
subdomains, all surrounded by spherical-shell subdo- 
mains that extend to the outer boundary of the com- 
putational domain (located far from the holes). The 
intermediate-zone cylindrical-shell and/or rectangular- 
block subdomains must overlap both the inner and outer 
spherical-shell subdomains in these grid structures to 
cover the computational domain completely. These grid 
structures are quite efficient, and have allowed our group 
to perform long stable inspiral calculations for a variety 
of simple binary black-hole systems [St], E^, [5^, [6(| , and 
also simulations of the merger and ringdown phases of a 
few of these simple cases p3L [56j . The overlap regions in 
these grid structures are well behaved in these successful 
cases. But in more generic inspiral and merger simula- 
tions, these overlap regions become significant sources of 
numerical error and instability. 

The overlap-related instability discussed above is most 
likely caused by the method we use to exchange informa- 
tion between adjacent subdomains. We set the incoming 




FIG. 1: Illustrates a basic element of the grid structure used 
to fill the 3D volume between a sphere and a concentric cube. 



characteristic fields, whose values are determined by in- 
terpolation from the adjacent subdomain, using a penalty 
method [6l|, [H, [H, [13] ■ This penalty method for impos- 
ing boundary conditions was derived explicitly for the 
case where adjacent subdomain boundaries touch but do 
not overlap. One possibility for resolving this instability 
would be to re-derive the appropriate penalty boundary 
terms for the case of overlapping subdomains. We have 
chosen the simpler option of constructing new grid struc- 
tures without overlapping subdomains. We do this by 
changing the intermediate-zone grid structure from the 
cylindrical-shells and/or rectangular blocks used previ- 
ously into a series of deformed rectangular blocks whose 
boundaries conform to the spherical shells used near the 
inner and outer boundaries. The new type of grid ele- 
ment needed for this construction is illustrated in Fig. [TJ 
which shows how the 3D volume between a sphere and a 
concentric cube can be mapped into six deformed cubes. 
This construction is based on the "cubed-sphere" coor- 
dinate representations of the sphere (i.e. the mapping of 
the six faces of the cube onto the sphere as illustrated in 
Fig. H]), that is widely used in atmospheric and geophys- 
ical modeling (6f| [6g, [67|, [68j]. Three-dimensional grid 
structures based on cubed spheres have also been used 
in other types of simulations, including a few black-hole 
simulations [69l. |70|. 

Our new grid structure is illustrated in Figs. [2] and [3] 
It consists of a series of spherical-shell subdomains and 
several layers of cubcd-sphcrc subdomains surrounding 
each black hole, as illustrated in Fig. [5] The two cu- 
bic blocks containing the black holes are surrounded by 
a series of rectangular-block subdomains to fill out the 
volume of a large cube. The center of this large cube 
is placed at the center-of-mass of the binary system, as 
illustrated in Fig.[3]for a system having a 2:1 mass ratio. 
This large cube is connected to the outer spherical-shell 



FIG. 2: Grid structure around each black hole. A series of 
spherical shells (two shells shown) is surrounded by a series of 
cubed-sphere layers (three shown) to fill up a cube. This fig- 
ure illustrates one quarter of the cube structure surrounding 
one of the black holes. 
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FIG. 3: Illustrates the intermediate grid structure surround- 
ing the two black holes, consisting of a set of rectangular 
blocks placed around the cubic block containing each hole. 
This collection of blocks filles up a cube which is surrounded 
by a series of cubed-sphere layers (two shown) all surrounded 
by a series of spherical shells (one shell shown) that are cen- 
tered on the center of mass of the binary system. 



subdomains by several additional layers of cubed-sphere 
subdomains, also illustrated in Fig.[3J 

Previous attempts to use this type of cubed-sphere grid 
structure in SpEC have always proven unsuccessful, due 
to numerical instabilities at the interdomain boundaries. 
These previous attempts used no spectral filtering, or a 
two-thirds anti-aliasing filter, in the cubed-sphere subdo- 
mains. It has been shown, however, that an appropriate 
filter is required for stability (and improved accuracy) of 
Chebyshev polynomial spectral expansions, such as those 
used in these cubcd-sphcrc subdomains. Such a filter 
must satisfy two important conditions. The first criterion 
is that when represented in physical space, this filter must 
act like a dissipation term in the evolution equations that 
vanishes on the subdomain boundaries (to avoid conflict- 
ing with the boundary conditions). The second criterion 
is that it must set the highest spectral coefficient to zero. 
Adding this type of spectral filtering turns out to be one 
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of the key elements in improving our numerical method 
enough to make possible the binary black-hole simula- 
tions described in Sec. Mil The needed spectral filter is 
applied to each field u(x, t) that is expanded as a sum of 
Chebyshev polynomials: 

t) = ^u k {t)T k {x). (10) 

k 

These are used for the radial-coordinate expansions in 
the spherical-shell subdomains, and for all three spectral- 
coordinate expansions in the rectangular-block and the 
cubed-sphere subdomains. Each spectral coefficient Uk in 
these expansions is filtered according to the expression, 

T{u k ) =u k er a{k ' k ^' 2 \ (11) 

after each time step. This filter satisfies the first nec- 
essary filter criterion by adding a 2p th -order dissipation 
term to the evolution equations which vanishes on each 
subdomain boundary [64[. For the binary black- hole 
simulations presented here, we use the filter parameters 
a = 36 and p — 32. The choice a = 36 guarantees 
that this filter satisfies the second necessary filter crite- 
rion by ensuring that the highest spectral coefficient is 
set to (double precision) zero, while the choice p = 32 
insures that this is a very mild filter that only influences 
the largest few spectral coefficients. This new filter is ap- 
plied wherever Chebyshev expansions, Eq. (|10j) . are used. 
We use the same filter we have used in previous binary 
black-hole simulations in the angular directions in the 
spherical-shell subdomains: we set to zero the top four £- 
coefficients in the tensor spherical- harmonic represcntion 
of each dynamical field [43, H|| . 

C. Adaptive Conforming Grid 

One of the important developments that allowed our 
group (several years ago) to begin performing successful 
binary black-hole inspiral simulations was the introduc- 
tion of the dual-coordinate-frame evolution method [47j . 
This method uses two distinct coordinate frames: One is 
a non-rotating and asymptotically Cartesian coordinate 
system used to construct the tensor basis for the compo- 
nents of the various dynamical fields evolved by our code. 
The second is a coordinate system chosen to follow (ap- 
proximately) the motions of the black holes. We fix the 
computational grid to this second coordinate frame, and 
solve the evolution equations for the "inertial-frame" ten- 
sor field components as functions of these "grid-frame" 
coordinates. The map M. connecting these grid frame 
coordinates x l to the inertial frame x l , can be written as 
the composition of more elementary maps: 

M=M K °M S . (12) 

A4k represents a "kinematical" map that keeps the cen- 
ters of the black holes located (approximately) at the 



centers of the excised holes in our grid structures (cf. 
Sec. Ill B| ). and Ms is a "shape-control" map that makes 
the grid conform (approximately) to the shapes of the 
black holes. The third major technical development 
(which makes it possible now for us to perform ro- 
bust merger simulations) consists of improvements in the 
choice of the shape-control map, .Ms, and improvements 
in the way this map is adapted to the dynamically chang- 
ing shapes of the black holes. 

The kinematical map A4k can itself be decomposed 
into more elementary maps: 

Mk=Mt°Me°Mr. (13) 

The maps AIt, -Me, and .Mr each move the centers 
of the black holes, but do not significantly distort their 
shapes: The map A4t translates the grid to account for 
the motion of the center of the system due to linear mo- 
mentum being exchanged with the near field and 
being emitted in gravitational radiation. The map .Me 
does a conformal rcscaling that keeps the coordinate dis- 
tance between the centers of the two black holes fixed in 
the grid frame, as they inspiral in the inertial frame. And 
.Mr rotates the frame so that the centers of the two black 
holes remain on the x axis in the grid frame, as they move 
along their orbits in the inertial frame. These maps are 
also used during inspiral simulations in the same way we 
use them for our merger simulations. So here we focus on 
the new shape-control map Afs, which is one of the crit- 
ical new developments that made the merger simulations 
reported in Sec. IIIII possible. 

The interiors of the black holes are excised in our sim- 
ulations at the spherical boundaries shown in Figs. [2] and 
[3] This is possible because, if this excision boundary is 
chosen wisely, the spacetime inside this boundary can 
not influence the spacetime region covered by our compu- 
tational grid. For hyperbolic evolution systems, like the 
generalized harmonic form of the Einstein equations used 
in our code [45j |. boundary conditions must be placed 
on each incoming characteristic field at each boundary 
point. Apparent horizons are surfaces that are often used 
to study black holes, because they can be found numeri- 
cally in a fairly straightforward way, and because (if they 
exist at all) they are always located within the true event 
horizons. If the excision bounarics for our computational 
domain were placed exactly on the apparent horizons, 
then all the characteristic fields of the Einstein system 
would be outflowing (with respect to the computational 
domain) since apparent horizons are spacelike (or null) 
hypersurfaces. Boundary conditions would not be needed 
on any field at such boundaries. Unfortunately we are not 
able to place excision boundaries precisely on the appar- 
ent horizons. So the best that can be done numerically is 
to place them slightly inside the apparent horizons (i.e., 
the apparent horizons must remain within the computa- 
tional domain) , if we are to avoid the need for boundary 
conditions on these excision surfaces. 

However, if an excision boundary is placed inside an 
apparent horizon, the outflow condition is no longer auto- 
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matic or simple: the condition depends on the shape and 
the location of the excision boundary, its motion with 
respect to the horizon, and the gauge. In our simula- 
tions, the excision boundaries are kept somewhat inside 
the apparent horizons, so the outflow condition can — and 
does — fail if we are not careful. One reason (and proba- 
bly the principal reason) we need to control the shapes 
of the computational domains through the map A^s is to 
keep pure outflow conditions on all the dynamical fields 
at the excision boundaries. We do this, as described in 
more detail below, by requiring that the excision bound- 
aries closely track the shapes and sizes of the horizons. 

Another (probably secondary) reason that shape and 
size control are needed in our numerical simulations is 
related to finite numerical resolution. If the excision 
boundaries had different shapes than the horizons, then 
some points on the boundaries would be located deeper 
into the black-hole inerior and hence closer to the space- 
time singularity. Higher numerical resolution would be 
needed at these points, and for fixed resolution, the 
amount of constraint violation and errors in the solution 
would be largest there. In many situations we find that 
numerical instabilities at these points cause our simula- 
tions to fail. This mode of failure can be eliminated by 
keeping the shapes and sizes of the boundaries close to 
those of the horizons. 

We have decomposed the map M, which connects 
the grid-frame coordinates x 1 with inertial-frame coor- 
dinates x l , into kincmatical and shape-control parts: 
M = Mk ° Ms- ^ will be convenient to have a name 
for the intermediate coordinate system whose existence is 
implied by this split. The map Mk connects the inertial- 
frame coordinates x % with coordinates x 1 in which the 
centers of the black holes are at rest (approximately). So 
it is natural to call these intermediate-frame coordinates, 
x\ the "rest frame" coordinate system. The shape con- 
trol map Ms connects the grid-frame coordinates x 1 with 
these rest- frame coordinates x l . 

It is useful to express the shape-control map Ms as 
the composition of maps acting on each black hole indi- 
vidually: Ms = Ms! ° Ms 2 - Such individual black-hole 
shape-control maps can be written quite generally in the 
form 

o A = e A , (w) 

4>A = 4>A, (15) 
fA = f A - f A {f A ,6 A ,4>A) 

fmax I 

X E E >?T{m m {eA,4>A), (i6) 

1=0 m=~l 

where (fA, 6~a,4>a) and (?a, 6a,<Pa) are, respectively, the 
grid-frame and rest-frame spherical polar coordinates 
centered at the (fixed) grid-coordinate location of black 
hole A = {1,2}. The functions Ja^a^a^a) are fixed 
functions of space in the grid frame. We note that these 
maps become the identity whenever fA(fA,0A,4>A) = 0, 
so making Ja^a, 9a, 4>a) vanish as f a —> oo ensures that 



the distortion is limited to the neighborhood of each black 
hole. The parameters X e ™(t) specify the angular struc- 
ture of the distortion map. They are determined by a 
feedback control system (discussed in Appendix p) that 
dynamically adjusts the shape and overall size of the grid 
frame coordinates relative to the rest-frame coordinates. 

The Caltech/ Cornell collaboration has used shape- 
control maps of this form, Eqs. fl4" ]) -(116 fl . in all of our 
recent binary black-hole simulations |53l . l54l |5(| . In 
those cases the functions Ja^a^a^a) were taken to 
be smooth functions of fA alone: roughly constant near 
each black hole and falling off to zero rapidly enough that 
A^Si and Ms 2 approximately commute. These functions 
had large gradients which made it difficult to control the 
grid distortion near the horizons without introducing ad- 
ditional unwanted distortions elsewhere. As a result, we 
were never able to achieve very robust shape control us- 
ing these maps. 

One of our major breakthroughs, leading to the suc- 
cessful merger results reported here, is an improvement 
in our choice of the map functions fA(fA,&A,4>A)- We 
now use simpler functions $a{ta, 6a, 4>a) that have much 
smaller gradients and that are exactly zero outside (non- 
intersecting) compact regions surrounding each black 
hole. This means the new maps A^Si and Ms 2 commute, 
exactly. These new maps can be defined everywhere be- 
tween the black holes because of the new non-overlapping 
grid structure discussed in Sec. Ill Bt but these improve- 
ments are achieved at the expense of smoothness: the 
new 6a, 4* a) arc smooth except at subdomain 

boundaries where they may only be continuous, not dif- 
ferentiable 1 . Fortunately, this lack of smoothness at 
subdomain boundaries is not a problem for our basic evo- 
lution method. In our multidomain code the equations 
are solved independently on each subdomain, and sub- 
domains communicate only by equating the appropriate 
characteristic fields at their mutual boundaries. These 
boundary conditions require no differentiation, so the 
grid coordinates themselves need not be smooth across 
these boundaries. 

The new ^a(ja, 6a,4>a) are chosen to be unity inside a 
sphere of radius oa centered on black hole A={i,2}, then 
decrease linearly with radius along a ray through the cen- 
ter of the black hole, and finally vanish on the surface of 
a centered cube of size 2a a- A two-dimensional sketch 
of this grid-frame domain is shown in Fig. 01 which is 
merely an abstract version of the type of grid structures 
we use around each black hole as illustrated in Figs. Q] 
and [2j This function Ja if a , 6 a > ) can be expressed 



1 The idea of introducing non-smooth time-dependent mappings 
was also suggested by Larry Kidder. 
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FIG. 4: Two-dimensional illustration of the grid-frame coor- 
dinate domain on which the function /a(7\4, Oa, 4>a) is defined 
in Eq. (TT7|). 

analytically as 

r A <b Al 
l A ~ P - A , b A <f A < p A , (17) 

OA - PA 

0, pa < r A , 

where 

Pa{0a,(/)a) = a A (18) 
x [max ( | sin #a cos U sin #a sin <^ | J cos #a | ) ] 1 

is the value of f A on the surface of the cube of size 2a a 
centered on black hole A={i,2}. The different cases of the 
maxima which appear in the denominator of Eq. Q18p 
correspond to the different faces of this cube. 

We note that the choice of Ja^a, 9a,4>a) in Eq. (fT?]) 
implies that the shape-control map Ms, an d hence the 
map between grid-frame and rest-frame coordinates is 
not smooth. Since the kinematical map Mk is smooth, 
this also imples that the full map M relating grid-frame 
to inertial-frame coordinates is not smooth either. This 
non-smoothness docs not cause a problem for our basic 
evolution method, but it implies that calculations that 
require smoothness must not be performed in the grid 
frame. For example, the representation of a smooth ap- 
parent horizon in grid coordinates will not be smooth. 
So apparent horizon finding and other calculations that 
require smoothness arc done in the smooth rest-frame or 
inertial-frame coordinates. 

The angular structures of the shape-control maps de- 
fined in Eq. (|14p - (|16p are determined by the parame- 
ters X 1 ™ (t) . These are chosen dynamically to ensure that 
the shapes of the excision boundaries match the evolving 
shapes of the apparent horizons. These maps also con- 
trol the sizes of the excision boundaries by adjusting X A °, 
which are chosn to ensure that the excision boundaries 
remain close to but safely inside the apparent horizons. 
FigurcOillustrates the effect of this distortion map on the 
structure of the grid surrounding a black hole in one of 




FIG. 5: Illustrates the inertial frame representation of the grid 
structure around one black hole, just at the time of merger. 
This grid structure has been distorted in relation to the iner- 
tial frame by the shape-control maps Ms described here. 

the generic merger simulations described in Sec. 11111 The 
choices of suitable target shapes and sizes for these shape- 
control maps, and the feedback control system that we 
use to implement these choices dynamically in our sim- 
ulations, are also among the critical technical advances 
that allow us now to perform robust merger and ringdown 
simulations. The details of exactly how this is done are 
somewhat complicated, however, so we defer their dis- 
cussion to the Appendix „ 

III. MERGER AND RINGDOWN 
SIMULATIONS 

In this section we present simulations of the merger and 
ringdown of binary black-hole systems using our Spectral 
Einstein Code (SpEC). In Section Ull Al we describe the 
binary black-hole initial data sets that we evolve. In Sec- 
tion IIII Bl we discuss our evolution algorithm, pointing 
out at which stage we employ the various improvements 
described in Section [TT] and what goes wrong when these 
improvements are not used. Finally, in Section IIII CI 
we describe six different binary black-hole merger and 
ringdown simulations peformed using these new meth- 
ods. We present snapshots of the shapes and locations 
of the horizons at various times during the merger and 
ringdown, and a demonstration of convergence of the con- 
straint violations. 



A. Initial Data 

Our initial data are in quasi-equilibrium [7ll . l72l . l73j 
(see also [74| [7f|), and are built using the conformal 
thin sandwich formalism ffdl , [77| with the simplifying 
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Case 


Mi /Mi 


Si /Ml 


§2 /Ml 


iVorbits 


A 


1 








16 


B 


2 








15 


C 


1 


-0.41 


-0.42 


11 


D 


1 


0.42 


0.42 


15 


E 


2 


0.2(£ -x)/y/2 


-0.4(1 + y)/V2 


8.5 


F 


2 


0.2{z-x)/V2 


-0A(z + y)/V2 


1.5 



TABLE I: Black hole binary configurations run through 
merger and ringdown using the methods presented here. For 
the initial spin parameters 5*1 /Mi and 5*2 /M$ , the z direction 
is parallel to the orbital angular momentum. 



choices of conformal flatness and maximal slicing. Quasi- 
cquilibrium boundary conditions are imposed on spheri- 
cal excision boundaries for each black hole, with the lapse 
boundary condition given by Eq. (33a) of Ref . [73[ . The 
spins of the black holes are determined by boundary con- 
ditions on the shift vector at each excision surface [72| . 

This formalism for constructing initial data also re- 
quires the initial radial velocity v r of each black hole to- 
ward its partner and an initial orbital angular frequency 
fio (chosen to be about the z axis without loss of gen- 
erality). These parameters determine the orbital eccen- 
tricity of the binary, and can be tuned by an iterative 
method [56|, [6(| to produce data with very low eccentric- 
ity. This has been done for most of the initial data sets 
described here (cases A through D in Tabic HJ), but is un- 
necessary for the purposes of the present paper, since our 
goal here is simply to document our improved methods 
for evolving binaries through merger. For cases E and F, 
v r and Qq are chosen roughly by using post-Newtonian 
formulae that ignore spins. 

Table U shows the mass ratio and initial spins for the 
black hole binary configurations we evolve here. The ini- 
tial data and first 15 orbits of inspiral for case A are 
identical to those presented in [49|, [53| . The initial data 
and first 9 orbits of inspiral for case C are identical to 
those presented in (56|. Cases E and F are fully generic, 
with unequal masses and unaligned spins, and exhibit 
precession and radiation-reaction recoil. 

B. Evolution procedure 

In this section we summarize our procedure for evolv- 
ing black hole binary systems through merger and ring- 
down, concentrating on the improvements discussed in 
Section [Til This procedure can be divided into several 
stages, beginning with the early inspiral and extending 
through ringdown: 

1. Evolve early inspiral using "quasi-equilibrium" 
gauge. 

2. Transition smoothly to damped harmonic gauge. 

3. Eliminate overlapping subdomains. 



4. Turn on shape control. 

5. Replace remaining inner spherical shells with cubed 
spheres. 

6. Turn on size control just before holes merge. 

7. After merger, interpolate to a single-hole grid to 
run through ringdown. 

The above ordering of these stages is not the only possi- 
ble choice: we have exchanged the order of some of these 
stages without trouble. For example, for the runs de- 
scribed here we perform stage [3] at the instant we start 
stage [5] The transitions between these stages are rel- 
atively simple, require little or no fine tuning, and can 
be automated. We now discuss each of these stages in 
turn. For several of these stages, we will illustrate the ef- 
fects of new improvements (Section [TTJ) for one particular 
evolution, case F of Table [TJ 

1. Early inspiral 

The gauge is chosen early in the evolution following the 
procedure of Ref. [4!| : the initial data are constructed in 
quasi-equilibrium, so we choose the initial H a to make the 
time derivatives of the lapse and shift zero. We attempt 
to maintain this quasi-equilibrium condition by demand- 
ing that H a remains constant in time in the grid frame 
in the following sense: we require 

d- t H- a = 0, (19) 

where H a is a tensor (note that H a is not a tensor) 
defined so that H a = H a in the inertial frame. The 
bars in Eq. (fH)]) refer to the grid frame coordinates. 
We refer to this as "quasi-equilibrium" gauge, although 
strictly speaking this gauge maintains quasi-equilibrium 
only when the spin directions are constant in the grid 
frame. 

For cases A-D, the simulations follow many orbits us- 
ing this quasi-equilibrium gauge. For cases E and F, in 
which the spins change direction in the grid frame, this 
gauge is not appropriate and causes difficulties, so for 
these cases we transition to harmonic gauge H a = very 
early in the inspiral. 

During inspiral, the map Me uniformly contracts the 
grid so that the grid-coordinate distance between the 
centers of the two black holes remains constant. Be- 
cause this contraction is uniform, the spherical-shell sub- 
domains inside each apparent horizon shrink relative to 
the horizon. This means that the excision boundary (the 
inner boundary of the innermost spherical shell) inside 
each black hole moves further into the strong-field region 
in the interior, towards the singularity. If this motion 
continues unchecked, gradients on the innermost spheri- 
cal shell increase until the solution is no longer resolved. 
However, the inner boundary of the next-to-innermost 
spherical shell is also shrinking, and it eventually shrinks 
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FIG. 6: The minimum of the outgoing (into the hole) charac- 
teristic speeds on the nine innermost spherical grid boundaries 
(labeled Ao . . . As) around the larger black hole, for part of 
run F. 



enough that the characteristic fields on this boundary 
all become outflowing (into the black hole); when this 
happens, the innermost spherical shell no longer influ- 
ences the exterior solution, so we simply drop that shell 
from the domain. This shell-dropping process occurs au- 
tomatically as successive spherical shells shrink relative 
to the horizons, This process is illustrated in Figure [51 
which shows the minimum of the (into the hole) charac- 
teristic speeds on spherical grid boundaries around one 
of the black holes. At any given time, the innermost of 
these boundaries is an excision boundary. All character- 
istic speeds at all points on the excision boundary (and 
hence the minimum) must be positive for the excision al- 
gorithm to be well-posed. The innermost shell is dropped 
whenever the characteristic speeds are positive on the 
inner boundary of the next shell. For example, surface 
A3 in Figure [6] is the excision boundary at t — 2AM, 
but after the characteristic speeds on A4 become posi- 
tive (around t = 26M), the innermost shell is dropped 
and A4 becomes the excision boundary. Similarly, A§ 
becomes the excision boundary at t = 38, and so on. It 
is important that the characteristic speed on surface A n 
in Figure [5] does not become negative before the speed 
on surface A n+ \ becomes positive; otherwise this shell- 
dropping procedure will fail. Shell dropping works well 
during most of the evolution but it fails in the last stage 
before merger, which is discussed in Section [III B 61 

When describing mergers, it is useful to introduce a 
measure of distance between the black holes. One such 
measure of distance is indicated on the horizontal axis at 
the top of Figured! This is the spatial (in each time slice) 
proper separation between the apparent horizon surfaces, 
obtained by integrating along the rest-frame x axis (recall 



that in the rest frame, the centers of the horizons always 
remain along this axis). This is not the true proper sepa- 
ration between the horizons because we do not minimize 
over all possible paths between all possible pairs of points 
on the two surfaces; however, we use this only as an ap- 
proximate measure of how far the black holes are from 
merger. Note that the horizons touch each other when 
(or possibly before) this proper separation measure falls 
to zero. 

It turns out that the improvements discussed in Sec- 
tion |TT] arc unnecessary during early inspiral for cases 
A-D and are not yet used at this stage. For instance, 
the domain decomposition uses overlapping cylinders and 
spherical shells, and not the cubed-sphere subdomains 
described in Sec. Ill Bl Wc find no problems while the 
holes are far enough apart that they remain roughly in 
equilibrium. However, the binary eventually becomes ex- 
tremely dynamical, and the black holes become signifi- 
cantly distorted. Unless the algorithm is modified, the 
simulation fails (typically because of large constraint vi- 
olations) before the black holes merge. 

Typically, it is possible to extend the quasi-equilibrium 
inspiral stage until the proper separation decreases to 
about 7M or 8M, but the next stages, Sees. IIHB 31 
and IIII B 21 can be done sooner if desired. (For exam- 
ple, cases E and F use a non-overlapping cubed-sphere 
grid from t — 0.) For the simulations presented here, the 
length of this quasi-equilibrium inspiral stage is variable. 
For instance, for case A in Table [T] this stage lasts for 15 
orbits, but for case F we transition to damped harmonic 
gauge and eliminate overlapping subdomains starting at 
t = because the black holes are initially very close to- 
gether. 



2. Transition to damped harmonic gauge 

When the inspirling black holes reach a proper sepa- 
ration of about 8M, we smoothly transition from quasi- 
equilibrium gauge to damped harmonic gauge at t = t g 
(g stands for "gauge") by choosing 



H a (t) = H a (t) 



+VL log Nfj t a - (IsN-^aiN 1 



(20) 



where H a (t) is the value of H a (t) obtained from the 
quasiequilibrium gauge condition (fT9)) . and a g is a con- 
stant. The last line in Eq. (|20[) is the damped harmonic 
gauge condition, Eq. (|5|). The values of [ij, and /ig are 
set according to Eq. ©, using 



fJ-o 



0, 



t < 



1 



-(t-t d y/*i } t>t d , 



(21) 



where td and a d are constants. It is necessary to choose 
<7g and ad sufficiently large or else the gauge becomes 
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unnecessarily dynamical, possibly resulting in incoming 
characteristic fields on the excision boundary. 

For all runs shown here, td corresponds approximately 
to a proper separation of 8M. We use the values a g = 
20M, ad = 50 M for the equal mass binaries (cases A, C, 
and D), a g = 15M, ad = 100M for case B, and a g = 
25M, ad = 50M for case F. For these cases, t g = td- For 
case E, the gauge is rolled off to harmonic early in the 
inspiral, so at t = td we need only to roll on the damped 
harmonic gauge, and we use ad = 40M. The constants 
in Eqs. (|20|) and (j2Tj) can be chosen quite flexibly: the 
only constraint is that a g and ad must be large enough 
to avoid significant spurious gauge dynamics, but small 
enough so that the simulation is using damped harmonic 
gauge before the black holes approach each other too 
closely. The a g and ad can always be made longer by 
choosing earlier transition times t g and td- 



3. Eliminate overlapping subdomains 

Although the domain decomposition consisting of over- 
lapping spherical shells, cylinders, and cylindrical shells 
is quite efficient for the early inspiral, it often suffers from 
numerical instabilities, particularly as the black holes ap- 
proach merger. Therefore, we regrid onto the cubed- 
spherc subdomains described in Sec. Ill Bl Use of these 
non-overlapping subdomains eliminates these remaining 
instabilities, and effectively increases our resolution. For 
most of the runs described here, we happen to regrid at 
the same time td that we begin the transition to damped 
harmonic gauge. However, for runs E and F, the entire 
simulation uses the non-overlapping cube-sphere subdo- 
mains. As far as we know, there is no reason that regrid- 
ding cannot occur at times other than td or t g . 



4- Shape control 



Proper separation/M 
7 6 5 



The next stage is to turn on shape control by introduc- 
ing the map A^s defined in Eqs. (fT^lTTBj) . The parameters 
A^™ that appear in this map are determined (via a feed- 
back control system) by Eqs. (|A.3|) and (|A.6|) for i > 0, 
and A™ is set to zero. This map deforms the grid so that 
the boundaries of the spherical shells inside the horizons, 
including the excision boundaries, are mapped to closely 
match the horizons' shapes. 

We do not turn on shape control until time t g + 2.5a g , 
which is when the coefficient multiplying the gauge term 
H a (t) in Eq. (|20[) becomes smaller than double-precision 
roundoff. This ensures that H a remains a smooth func- 
tion of space in inertial coordinates. If shape control 
were turned on too early, H a would fail to be everywhere 
smooth because the quasi-equilibrium-gauge quantity H a 
is a smooth function of space in the grid-frame coordi- 
nates (and its numerical value is constant in time for 
each grid point), and the deformation is only C° across 
subdomain boundaries. 




FIG. 7: The minimum of the outgoing (into the hole) char- 
acteristic speeds on the remaining four spherical-shell bound- 
aries (labeled A5 . . . As) around the larger black hole, for a 
portion of run F. The dashed curves correspond to runs with- 
out shape control, and the solid curves correspond to runs 
with shape control turned on at t = 62M. The run with- 
out shape control terminates at proper separation 5.2M (time 
t = 106M) because the excision algorithm fails on surface A7. 



Shape control is necessary for the shell-dropping pro- 
cedure described in Scction llll B ll to remain successful as 
the holes approach each other. For a shell to be dropped, 
the entire inner boundary of that shell must remain an 
outflow surface until the entire inner boundary of the 
next shell becomes an outflow surface. If these bound- 
aries do not have the same shape as the horizon, then 
typically some portion of either boundary moves too close 
or too far from the horizon, and violates this condition. 
An example is shown in Figure [7] Like Figure [SJ Figure [7J 
shows the minimum characteristic speeds on the bound- 
aries of inner spherical shells for a portion of case F, but 
Figure [Jj shows two separate evolutions: one with shape 
control turned on at t — 62M (solid curves) and one 
with no shape control (dashed curves). Without shape 
control, the minimum characteristic speed on surface A-j 
stops growing around t — 100AZ and becomes negative 
soon thereafter. Because Aj is the excision surface at 
t ~ 100A/, the excision algorithm becomes ill-posed and 
the simulation stops at t = 106M, well before merger 
(which occurs at t ~ 128M). In the simulation with ac- 
tive shape control, the same boundary Aj shows no such 
problems — in fact it is dropped at t = 102M when all 
speeds on the next shell, Ag become positive. 
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FIG. 8: Constraint norm ||Cgh|| as a function of time for a 
portion of case F. The dashed curve corresponds to a run with 
spherical shells around each hole, and the solid curve corre- 
sponds to an identical run with these inner spherical shells 
replaced by cubed-sphere shells at t — 80M. The constraint 
norm grows larger and begins to grow earlier when spherical 
shells are used. 



5. Eliminate inner spherical shells 

Spherical-shell subdomains are extremely efficient 
when the solution is nearly spherically symmetric, be- 
cause these subdomains use Yi m basis functions, which 
are well-suited for nearly spherical functions. However, 
once the holes become sufficiently distorted, it is difficult 
to resolve the region near each horizon using spherical 
shells. At this point in the evolution, we replace the in- 
ner spherical shells with cubed-sphere shells. Figure [8] 
shows the constraints as a function of time for two evo- 
lutions of case F: one with inner spherical shells (dashed 
curves), and another in which the inner spherical shells 
were replaced with cubed-sphere shells at t — 80M (solid 
curves). The quantity plotted is ||Cc7ff||, the L 2 norm 
over all the constraint fields of our first-order generalized 
harmonic system, normalized by the L 2 norm over the 
spatial gradients of all dynamical fields (see Eq. (71) of 
Ref. [H|). The L 2 norms are taken over the portion of 
the computational volume that lies outside the apparent 
horizons. For the spherical-shell case in Figured the con- 
straints in the spherical shells grow large enough that the 
simulation terminates at t = 129. 5AT, when the proper 
separation has fallen to 0.5AT. Note that these simu- 
lations typically proceed through merger even if we do 
not replace spherical shells with cubed-sphere shells. So 
while this replacement step is not strictly necessary, elim- 
inating spherical shells improves the accuracy of the sim- 
ulation (as measured by the constraint quantity ||CGff||) 
by almost an order of magnitude. 



6. Size control 

Eventually, even using shape control, shell dropping 
typically fails before the horizons merge: the outflow con- 
dition on the inner boundary of the inner shell fails before 
the outflow condition on the inner boundary of the next 
shell becomes satisfied. This problem occurs because the 
map Me shrinks the entire grid — including regions in- 
side each hole — faster and faster as the holes approach 
each other. As soon as the velocity of excision bound- 
ary towards the center of the hole becomes large enough, 
the boundary becomes timelikc, and the excision algo- 
rithm fails. To remedy this, we turn off shell dropping 
and we turn on size control. Size control slows down 
the infall of the excision surface (by pulling the excision 
sphere towards the horizon), thus keeping the character- 
istic speeds from changing sign. Size control is imple- 
mented by changing the map parameters X A °(t), which 
previously were set to zero, to values given by Eqs. (|A.3|) 
and (jA.lOp . This causes the size of each excision bound- 
ary to be driven towards some fraction rj of the size of 
the appropriate horizon. For the equal-mass runs A, C, 
and D we use 77 = 0.8 for both black holes. For the 
unequal-mass run B and the generic cases E and F we 
use 77 = 0.9 and ?y = 0.7 for the smaller and larger black 
holes, respectively. These values do not require fine- 
tuning; for instance, case E still runs through merger and 
ringdown if we change 77 to 0.9 for both black holes. The 
main criterion for choosing this parameter is that a grid 
sphere of radius b~A (see Figure 0]) should not be mapped 
outside the cube of side 2a a- If this occurs, the map 
A^s becomes singular and the run fails. Figure O shows 
the minimum characteristic speed at particular domain 
boundaries for two evolutions of case F: one with size 
control and one without. Without size control, the out- 
flow condition on the excision boundary fails at proper 
separation 1.67A/ and time t = 127.24M. With size con- 
trol, the evolution proceeds to proper separation 0.015A/, 
time t = 130.17, well past the formation of a common ap- 
parent horizon, which occurs at proper separation 1AM 
(time t = 128.23AT ). 



7. Ringdown 

At some time t = t m shortly after a common appar- 
ent horizon forms, we define a new grid, composed only 
of spherical shells. This grid has only a single excision 
boundary inside the common horizon. We also define a 
new map A4 ringdown : x 1 — > x l that can be written as a 
composition of simpler maps: 

Af ring down : = Mrf ° Mei ° Mr o M S3 . (22) 

This is similar to the ringdown maps described in 
Refs. (53l . [Hi| (which had no translation) and Ref.psij 
(which had no rotation). The new translation map A^Tf 
is tied to the center of the new common horizon, and 
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FIG. 9: The minimum of the outgoing (into the hole) charac- 
teristic speeds on the two innermost cubed-sphere boundaries 
(labeled Bo and Bi ) around the smaller black hole, for a por- 
tion of run F. The dashed curves correspond to runs without 
size control, and the solid curves correspond to runs with size 
control turned on at t = 120. 2M. Without size control, the 
innermost cubed-sphere is dropped at t = 124. 5M and Bi 
becomes the excision boundary. However, the characteristic 
speeds become negative on Bi at t = 127M and excision 
fails. With size control, the characteristic speeds on the ex- 
cision boundary Bo remain positive until proper separation 
0.075M, well after merger. 



is not continuous with the old translation map Mt ex- 
cept near the outer boundary where both translation 
maps are the identity. The new expansion map A^Ef 
is chosen to be continuous with the old expansion map 
A4e at the outer boundary, but it is the identity near 
the merged black hole. This new expansion map A^e£ 
smoothly becomes constant in time shortly after merger. 
The rotation map Mb, is continuous at merger, and af- 
ter merger it smoothly becomes constant in time. The 
map A4s 3 has the same form as Eqs. (fTlHTTj) except the 
function pa(9a,4>a) defined in Eq. (jT5J) is replaced by 
f>3 = constant, where the constant value corresponds to 
a subdomain boundary. 

As described in Ref. [53|, at t = t m the parameters 
\ ™ are chosen so that the common apparent horizon 
is stationary and spherical in the grid frame. Similarly, 
the map Mti is chosen at t — t m so that the center 
of the horizon is located at the grid-frame origin and is 
stationary in the grid frame. 

Once we have defined the maps at t = t m , we then in- 
terpolate all variables from the old grid to the new grid. 
Note that the grid frame changes discontinuously in time 
at t = t rn . Because of this, we take care to properly 
use both the new map -M r i ng down and the old map M. in 
doing this interpolation, so that the inertial frame and 



FIG. 10: Three snapshots of the apparant horizons of the 
non-spinning equal- mass binary black hole merger (case A). 
Solid curves are the orbital plane cross section when a com- 
mon horizon is first detected; dotted curves represent the time 
of transition between the binary merger and the single hole 
ringdown evolutions; dashed curve shows the final equilibrium 
horizon. 

quantities defined in that frame remain smooth. In par- 
ticular, the gauge source function H a is still determined 
by Eq. (|20p during ringdown, and remains smooth at 

t tm . 

For t > t m , the parameters A^™(t) are determined by a 
feedback control system that keeps the common apparent 
horizon stationary in the grid frame. Likewise, the map 
MTf keeps the common apparent horizon centered at the 
origin in the grid frame. 



C. Results 

In this section we present some results from the sim- 
ulations listed in Table |TJ We first show snapshots of 
the inertial-coordinate shapes of the apparent horizons 
during merger and after ringdown. Figures 1101 1111 1121 
and H3] show cross-sections of the horizons in the orbital 
plane for cases A-D; in these cases, the orbital plane is 
well-defined and is constant in time. For all cases, we 
show the apparent horizons of the individual black holes 
and the common apparent horizon at the time when the 
common horizon is first detected (solid curves), and at 
the time t m when we transition to a new grid that has a 
single excision boundary (dotted curves). We also show 
the apparent horizon of the final remnant black hole after 
it has reached equilibrium (dashed curves). Figures [Ml 
and Qj)] show cross-sections of the horizons for cases E and 
F, in which the orbital plane precesses. For these cases 
we show cross-sections in the coordinate plane (defined 
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FIG. 11: Apparent horizon snapshots as in Figure llOl except 
for a 2:1 mass ratio non-spinning binary black hole merger 
(case B). 




FIG. 12: Apparent horizon snapshots as in Figure llOl except 
for an anti-aligned-spin equal-mass binary black hole merger 
(case C). 

using the flat metric) that is perpendicular to the instan- 
taneous orbital angular velocity at the time the common 
horizon is first detected, and that passes through the co- 
ordinate center of the common horizon at this time. We 
show cross-sections at three different times: the time of 
horizon formation, the time t m , and a time 300M after 
merger, all with respect to the same plane. The rem- 
nant black hole in cases E and F have nonzero linear 
momentum both because of radiation reaction and be- 
cause of some nonzero linear momentum present in the 
initial data; this is why the centers of the horizons shown 
in Figures [HI and [Pol change with time. 

For case F, we also show the constraint norm as a func- 




FIG. 13: Apparent horizon snapshots as in Figure [TOl except 
for an aligned-spin equal-mass binary black hole merger (case 
D). 




FIG. 14: Three snapshots of the apparant horizons of a 
generic binary black-hole merger (case E). The plane of the 
cross-sections is a coordinate plane perpendicular to the in- 
stantenous orbital axis at the time the common horizon is first 
detected. The additional degree of freedom is fixed by having 
this plane go through the coordinate-center of the shape of 
the common horizon at its first detection. Solid curves are 
the cross section when a common horizon is first detected; 
dotted curves represent the time of transition between the bi- 
nary merger and the single hole ringdown evolutions; dashed 
curve shows the horizon at a time well after merger. 

tion of time in Figure 1161 We plot the same quantity 
||Cgh|| a s shown in Figure [H] This quantity is shown for 
four numerical resolutions, and is convergent at all times 
(although the convergence rate is smaller near merger 
when the solution is most dynamical). The constraints 
are largest at t = t m , when we transition to a grid with 
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FIG. 15: Apparent horizon snapshots as in Figure [Til except 
for the generic binary black- hole merger case F in Table [T] 
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FIG. 16: Constraint norm ||Cgh|| as a function of time for 
the generic (case F) binary black-hole merger and ringdown, 
computed with four different numerical resolutions. The small 
inset graph shows that numerical convergence is maintained 
(at a reduced rate) even during the most dynamical part of 
the merger. 



a single excision boundary. Just after t = t rn the con- 
straints decrease discontinuously by a small amount be- 
cause part of the computational domain has been newly 
excised. The approximate number of grid points for the 
four resolutions is {79 3 , 87 3 , 95 3 , 103 3 } just before merger 
and {44 3 , 51 3 , 57 3 , 63 3 } during ringdown. 



APPENDIX: CONTROLLING THE MAPS 

The purpose of the shape-control maps -Ms A , with 
A={i,2}, is to distort the grid structures around the black 
holes so the excision boundaries of the computational do- 
main are mapped to surfaces lying just inside and having 
the same basic shapes as the apparent horizons. Choos- 
ing the right maps is equivalent to choosing the right 
target values for the parameters A^ Ti (i) that define these 
maps, cf. Eq. Q16p . This Appendix describes in some 
detail how these parameters are chosen, and what target 
surfaces are used to fix these maps in the merger simula- 
tions described in Sec. IIIII 

The grid structures used in our merger simulations, cf. 
Sec. Ill Bl have excision boundaries that are grid-frame 
coordinate spheres. The target surfaces 5j to which we 
want to map these grid-frame spheres can be written in 
the form 



(Al) 



£=0 m=- 



where the expansion coefficients A 1 ™ (t) define the target 
surface and (ta^a^a) are rest-frame coordinates. 

Let f^ arget denote the radii of the grid-frame coordinate 
spheres that are to be mapped onto the target surfaces 
Sj. The rest-fr ame representations of these coordinate 
spheres, using Eqs. (Q3HH]), are: 



TA 



- target 



£=0 m=-- 



Note that the functions Ja that appear in Eq. (|T7j) are set 
equal to unity here because we always choose the target 
grid-spheres f^ algct to be smaller than the outer radii of 
the largest spherical subdomain layers: b a > f^ arset . It 
follows that the target spheres will be mapped to the 
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shapes of <5>T when 



A lm 
~ 1V A 



Y 00 
(*), 



= 0, 
> 0. 



(A.3) 



These conditions can not be imposed directly for a num- 
ber of reasons, but they can be enforced approximately 
using a feedback control system. The control system used 
in the merger calculations described in Sec. IIHI is the 
same one used to perform our earlier binary black-hole 
simulations [47j . so we will not describe it in detail here. 

To complete the specification of the shape-control 
maps we must choose the target grid-spheres f J; arset , and 



Prom Eq. (TO)) it fol- 



the target surface parameters A^™. 

lows that the choice of f^ argct only effects the t = size 
control part of the distortion map. So until size control 
is activated late in our runs, there is no need to specify 
what f^ arget actually is. We think of it as being the ra- 
dius of the suitably scaled apparent horizon, but this fact 
does not influence the shape control part of the map in 
any way. Once size control is activated late in the run, 
then we require f A algct to be the radius of the excision 
boundary: 



; target 



(A.4) 



Next consider the choice of target surfaces S A , start- 
ing with the I > contributions that control the shape 
but not the overall scaling of the map. The idea is to 
choose the target surfaces iSj to be similar in shape to 
the apparent horizons 'Ha- The H.a can be represented 
as smooth surfaces in the rest-frame coordinate system: 



1 for x 3> 1. The p-dcpcndcncc maintains these asymp- 
totic forms, but allows the transition between them to 
occur more quickly than the p = 1 case. 

We have found no need to introduce the I = size 
control parts these maps, until the final plunge phase 
just before the black holes merge. So during most of our 
merger simulations we simply set X A °(t) = 0. This al- 
lows the apparent horizons to grow relative to the grid 
coordinates as the map M. e contracts the rest-frame rel- 
ative to the inertial coordinates. This growth in the ap- 
parent horizons in the grid frame allows us to drop a 
number of spherical subdomain layers as the evolution 
progresses, and this helps remove unwanted constraint 
violations from the computational domain. 

At some point however, we may not be able to drop 
additional spherical subdomain layers: either because we 
run out of layers, or because the inner boundary ex- 
periences incoming characteristic fields when the appar- 
ent horizon expands too quickly during the final plunge. 
When either of these conditions occurs we turn on the 
i = part of the shape-control maps, to keep the ra- 
dius of the excision boundary f 
apparent horizon. It will be useful to define 



J| x close to the size of the 



Ar, 



vRa - fT- 



(A. 



We choose r\ so that when Ata < there is no need to 
impose size control, and use r\ in the range 0.7 < rj < 0.9 
for the merger simulations reported in Sec. IIIIl When 
Ata becomes positive we want to turn on size control, 
but we want to do it in a fairly smooth and continuous 
way. This is done by defining the function: 



ta 



E E 

£=0 m=-t 



S A m (t)Y e , 



(A.5) 



To keep the shapes of the target surfaces similar to the 
shapes of the apparent horizons, the target surface pa- 
rameters A^™ should be made proportional to S l A m . We 
find it is appropriate to scale these coefficients by a fac- 
tor, G(Ra), which depends on the average radius of the 
apparent horizon, Ra 



QOOy . 

& A Yqo- 



A e r(t) = G{R A )ST{t). 



(A.6) 



The larger the apparent horizon radius in relation to the 
desired excision radius, the smaller this scaling factor 
must be to maintain an appropriate shape for the ex- 
cision boundry. In practice, we find the scale factor 



G(R A ) 



OA_ 

Ra 



tanh p 




(A.7) 



works quite well for p = 2. This scaling factor is near 
unity when Ra < a a and decreases like 1/Ra for larger 
values of Ra ■ The tanh x function is introduced here be- 
cause it is linear for small values of \x\, and approaches 



P{Ar A ) = 



tanh 



0, 

/ A0Ar A 
V Ra 



Ar A < 
Ar A > 



(A.9) 



which vanishes for Ata < 0, and asymptotically ap- 
proaches P(Ar A ) — > 1 for Ar A > -Ra/40. We find that 
an appropriate value for the overall scale of the target 
surfaces S A , i.e. their i = components, to be 



A oo 



Ar A P(Ar A ) 



00 



(A.10) 



Once size control is activated, the target radius is given 



by Eq. (TOl) : 



; target 



This choice implies that 



the average rest-frame radius of the excision boundary, 
Rj? = ~ -^a^oo, approaches t)Ra when Ata > 
R A /40. 

To summarize, the appropriate shape and size con- 
trol of the inner grid structures in our binary black- 
hole merger simulations are deterined by the target 
shape-control surfaces S A , whose definitions are given in 
Eqs. (|A.6|) and (|A.10[) . The shape-control map needed to 
distort the grid coordinates to the shape of S A , is given in 
Eq. (|A.3[) . This mapping is enforced with an appropriate 
feedback control system. 
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